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Abstract 

Using a global equation of state, empirically derived by Luding, we ac- 
curately model the density profile of a two-dimensional hard sphere system 
with diameter D and mass m under gravity with a given temperature T 
[Physica A, 271, 192 (1999)]. We then compare our solutions to MD sim- 
ulated data. From the density profile, we can then solve for the critical 
temperature T c , which we define as the temperature at which the system 
begins to condensate. Then, if T is below T c , there is some frozen portion 
of the system. We derive a formula for the number of frozen layers fif, 
and compare our solution to the number of frozen layers in our simulated 
data. 

PACS number(s) 64.70.Dv, 05.20.Dd, 51.10.+y 

1 Introduction 

Because of the large number of relatively small particles that compose a 
granular system, granular media tends to act like a collection of microscopic 
molecules. However, individual grains in a granular system have a macroscopic 
mass. Thus, particle collisions and gravity play a much larger role than in molec- 
ular systems, making it impossible to fully explain their mechanics in terms of 
kinetic theory alone. Understanding hard sphere (HS) systems is crucial for an- 
alyzing important physical phenomena in our world, such as avalanches, earth- 
quakes, and powder production. Newer models that correct for characteristics 
such as higher densities and dissipation in granular materials have succeeded 
in explaining simulated and experimental results that were previously left un- 
solved. 

In a paper by Hong, [1], a method for modeling a granular system was devel- 
oped by observing the effect of lowering the excitation level of a highly excited 
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system of grains. The excitation of grains is modeled as a temperature reservoir 
in contact with the system. When the kinetic temperature, T = ^m(v 2 ), is 
high, particles in the system behave like a gas, where the mean free path is 
much larger than the particle size. The particles in this state are free to ex- 
change positions with neighboring particles. At the other extreme, T = 0, the 
particles remain motionless at the bottom of the container. In this state, the 
particles behave like a solid, remaining locked in a cage, not free to exchange 
positions with any neighboring particles. As the temperature is lowered from 
the gas-like state, we note that a portion of the particles begins to move to the 
bottom, forming a solid regime, where a layer of particles is locked in a cage. As 
the temperature gets closer to T = 0, the number of particles condensing to the 
bottom increases, causing the solid regime, and hence, the number of condensed 
layers, to grow. This observation lead to the derivation of a critical tempera- 
ture, T c , at which the solid regime begins to form. Hong was able to determine 
this critical temperature as a function of the control parameters of the system. 
He also showed that there is a way to predict the solid regime as a function of 
T once it has dropped below the critical temperature, T c . In another paper by 
Quinn et al [2], the theory proposed by Hong was tested and verified using an 
event driven(ED) molecular dynamics simulation code. The results of the ED 
simulations were compared with the proposed theory. Hong derived results in 
both two and three dimensions, choosing a particular correlation function for 
each. In two dimensions, Hong chose to use the correlation function proposed 
by Ree and Hoover [3], expressed as follows: 

1 - ai4> + a 2 (t> 2 
x{4>) = (l-acf,) 2 ' 

with a = 0.489351 • (tt/2) w 0.76867, a x = 0.196703 • (tt/2) « 0.30898, and 
a 2 = 0.006519 • (tt 2 /4) ~ 0.0168084. Using this correlation function, Hong 
derived the following expression, 

-P(C-fi) = ln0 + Cl ^ + c 2 ln(l-a(/ N ' 
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and 4>o is the density at £ = 0. Furthermore, Hong was able to determine the 
following expressions for the critical temperature T Cl 

_ mgD^p 

J-c — 

Mo 

and the number of frozen layers £/, 




The theory proposed by Hong broke the system of hard spheres into two 
regions, a solid region and a liquid region. The density profile proposed by 
Hong was shown by Quinn to only be valid for the liquid region of the system. 
This is because the correlation function proposed by Ree and Hoover, was not 
valid for a high density region of hard spheres. Therefore it makes sense that 
the density function derived from this correlation function, would be inadequate 
for modeling the high density of the solid region. In a paper by Luding et al, 
[4] a new correlation function was proposed for a two dimensional system. This 
correlation function is valid for both high and low density regions of a system 
of hard spheres. In this paper, we will derive a new density profile using the 
correlation function proposed by Luding, and compare it to the data produced 
by event driven molecular dynamics simulations. We will also derive expressions 
for the critical temperature, T c , and the number of frozen layers, Q. We will 
also compare these expressions to the simulated data. 



2 Liquid-Solid Transition for Hard-Sphere Sys- 
tems 

In his paper, Hong started with the Enskog Equation for hard spheres [5], 

— + v • Vf - mg— = J E , 

at ov z 

where the Enskog collisional operator, Je, is given by 

Je = D 2 J d 3 Vl J d 2 c(e ■ g)[/(r, v')/(r + De, vi) X (r + De/2) 

-/(r,v)/(r-De,vi) X (r-£>e/2)]. 

From this equation, he derived the density profile used to describe the liquid 
region of the hard sphere system. The same expression can be derived start- 
ing with the pressure difference equation, taught in most introductory physics 
classes. The pressure difference equation that describes the change of pressure 
P over height z, in terms of gravity g and the density p, 

dP m 
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As a simple example, we can show how this simple equation can be used to 
derive the Maxwell-Boltzmann Distribution. Starting with the Ideal Gas Law, 



PV = nkT, 

where k is Boltzmann's constant and T is the temperature of the system, one 
can get the following Ideal Gas Pressure in terms of the density, p: 

P = pkT 

Substituting this expression into Eq.(l) and solving for p, we obtain the following 
density profile, 

p = p ex V (^f). 

This is the Maxwell-Boltzmann Distribution, used to model a gaseous phase in 
statistical mechanics. Since the function only applies for the gaseous phase, it 
does not accurately fit the liquid-solid transition in our system of hard spheres. 

From the kinetic theory, we can relate pressure and density, where m equals 
the mass of a single particle, 

P=^[p + 2p X (p)]. 

x(p) is a correlation function that takes into account the probability of a collision 
occurring among particles. 

In 2001, Stephan Luding [4] proposed the following Global Pressure Equation 
to combine the low-density pressure and high-density expression: 

P = p + pP 4 + pm(p) [P dense - Pi] (2) 

where 

Pi == -lp + 2pxi(p)l 

1 — — pTS 

Xi{p) = (T^# - 8tFpf> (3) 

and 

C 

Pdense = [1 + C\{p max — p) + C 3 (p max —/£>)]— 1 

Pmax P 

X4 is a correlation function, similar to the form proposed by Ree and Hoover. 
The first term in Eq.(3) is a simpler correlation function, introduced by Hender- 
son [6,7]. The correlation function is determined via a virial expansion around 
low densities. The value of P4, taken at contact, accounts for the excluded vol- 
ume effect and the increase of the collision rate with density. Pdense is Luding's 
proposed corrected high-density pressure with the numerically fitted constants 
Co w 1.8137, Ci = -.04, and C 3 = 3.25. 
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Included in the Global Pressure Equation, is the following empirical merging 
function, 

m(p) = 



1 + exp[ 



with the numerical fitting constant m = 0.0111. Unfortunately, this function 
inhibits our ability to analytically solve for the density function p. We can 
obtain an integral expression for the density profile by substituting Eq.(2) into 
the pressure difference equation, Eq.(l), and solving for the position z. We 
obtain: 



—T r p 1 

mg J Po p 



2p { C m(p) + C dm(p) + C m(p) | 

(pmax ~ P) 2 {Pmax - P) dz p{p m ax ~ P) 



/ x / x C C 3 dm(p) C C 3 m(p) 
2C C 3 (p max - p)m{p) + - / ' + , 2 + 

{Pmax - P) dz p(p m ax - pY 



{CoCi - l)ra(p) dP 4 dP 4 dm(p) 

~p 2f>m{p) ^z- 4m(/3) ^7 4pPi ^z~ dp] + Z ° (5) 

Unfortunately, there is no analytical solution to Eq.(5), but using numerical 
methods, one can obtain the solution, and hence, a single function to completely 
describe the density profile of a two dimensional system of hard spheres under 
gravity. 



3 The Critical Temperature and Number of Frozen 
Layers 

One can define the critical temperature as the point at which the density at the 
bottom layer, p D , becomes the closed packed density p c such that 

Po(T c ) = p c - 

For T < T c , a portion of the particles near the bottom settle into their minimum 
energy positions. The particles in this region are essentially locked in a small 
cage, where they are free to slightly wiggle, but are never able to exchange 
positions with other particles in the system. Hence, they form a crystal-like 
structure. 

The maximum density p a depends on the underlying crystalline structure. 
The density of square lattice packing is p — j , and the density of a hexagonal 

lattice is p = ZL g^. We say that particles become locked in position when the 
density reaches the square-packed density p c — j . However, we will see that the 
lattice of the solid regime can fluctuate between square and hexagonal packing. 

We can use Eq.(5) and particle conservation to derive an expression for the 
temperature, T c at which the solid region begins to form at the bottom. We 
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begin with the conservation of particles in a column which can be expressed as 
follows: 



/ z(p)dp = h, 

J Pn 



(i 

'Pa 

where h is the height of the column of balls at a particular temperature, T. The 
height h can be rewritten as the number of layers simply by factoring out the 
diameter of a particle, D, such that 



_ h 

This allows us to rewrite Eq.(6) as 



i r° 

- J z{p)dp. (7) 



'p. 

The variable (3 can be defined as 

a = mg 

P rp 

We can now define the constant 

= mg 

where T c is the critical temperature and the constant h Q such that 

h = Dp = ^. 

1 c 

Now Eq.(5) can be rewritten as 

-p c z = I(p) - I(p ), (8) 

where the function I(p) comes from the upper limit of the integrals in Eq.(5), 
and I{p ) comes from the lower limit of the integrals in Eq.(5). Then we can 
rewrite Eq.(8) as 



where f3h — I(p )- 
Solving Eq.(9) for z yields 



-0 c z - I(p) - 0h, (9) 



z = -j[I(fi)-p c h]- (10) 



We can now substitute Eq.(10) into Eq.(7) to obtain 



,= ^j\p)dp 
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PcD [ 





[I(p) - c h]dp 



[I(p)dp - c hp] 

Po 



= ito Jo ° [1{p)df ' Pclp] - 
However, by our previous definitions, we can state that 

h Q T c h D 

p= m = D^ 9 (12) 

Using Eq.(12), we obtain an expression for the critical temperature T c , 

T ° ~ — ' (13) 
where h is obtained using the following expression: 

ho= [ P °[I(p)-I(po)]dp. 
Jo 

Using the Global Equation of State and Eq.(5), we calculate the value of h a 
to be 26.8097. Thus, 

T c = ^i. (14) 
26.8097 v ' 

We can now derive an expression for the number of frozen layers, using the 
critical temperature from Eq.(14). The number of frozen layers, £/ will simply 
be the total number of layers, /i, less the number of fluidized layers on the top, 
h /D(i. This gives us the following expression for the number of solid layers: 

h a pD/3 c 

c ' = M -^ = ^-^r (15) 

Finally we obtain the following expression for the number of layers by sim- 
plifying Eq.(15): 

Cf=p(1-yJ- (16) 

The expression we have derived is exactly the same as that derived by Hong 
in [1], with a modified value of T c , due to the Global Equation of State. 
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4 Molecular Dynamics Simulations 



We have used an event driven(ED) molecular dynamics code to test the conden- 
sation of hard spheres due to gravity. We refer the readers to references [8,9,10] 
for details of the algorithm regarding the collisional dynamics. This particular 
code takes into account the rotation of the hard spheres as well as a way to han- 
dle the inelastic collapse. The difference between ED and soft-sphere molecular 
dynamics codes lies in the way time is advanced during the simulations. Instead 
of using finite time increments to advance the system, the ED algorithm finds 
and advances the system to the next possible event, usually a collision between 
particles or between a particle and a wall. This process advances the system in 
time, but by different time steps for each event that occurs. The thermal reser- 
voir of our system was modeled using a white noise driving, first introduced 
by Williams and Mackintosh [11]. In this model, a random velocity kick, Av, 
is added to each particle's velocity every set interval of time, At. These ran- 
dom velocity kicks are not correlated with each other in any way. The program 
allows us to set the temperature of the system by setting an input parameter 
controlling the width of the Gaussian from which the random kicks are drawn. 
In effect, the temperature parameter controls the temperature of the reservoir 
and hence, the kinetic temperature of the system. Note that we are not driving 
the system with a bottom wall connected to a temperature reservoir, which is 
used fairly often as a model for the vibrating bed. 

Fig.(l) displays data from a typical simulation of 1000 particles, with initial 
layers , fi — 40 layers, radii r — .0001 m, and mass m = 8.378 x 10~ 9 kg. 
This data is fit with the three density functions derived previously, the Mazwell- 
Boltzmann Distribtuion, the Enskog approximation using P 4 , and Eq.(5) derived 
from the Global Equation of State. It is clear that the Maxwell-Boltzmann 
Distribution only works for the very low density region, or the gaseous region as 
expected, and the Enskog approximation breaks down once the particles reach 
their maximum denisty, or solid-like region. Eq.(5) fits the data for all densities, 
making the Global Equation of State the most accurate for representing our 
granular system. 

We now present typical configurations as a function of temperature using 
Eq.(5). In Figs. (2), (3) and (4) there are three representative samples at different 
temperatures for 1000 particles with initial layer numbers of [i = 40, of mass 
m = 8.378 x 10~ 9 kg and diameter D = 0.0002 m. The initial layer thickness 
is [i = 40, and the gravitational constant g = 9.81 m/s 2 . Seen in Fig. (2), at 
a high temperature (Ti > T c and p < p c ), all the particles are fiuidized and 
dynamically active. The density profile is fit using the density derived from the 
Global Pressure Equation, Eq.(5). Since the condensation picture is based on 
elastic hard spheres, we used both a completely elastic system, where one burst 
of energy is used to start the system, as well as coefficients of restitution that 
were close to unity in these simulations. We note, however, that it was shown 
in [4] that the presence of dissipation does not change the condensation picture. 
The critical temperature is just shifted linearly on the order of the coefficient 
of restitution. In our simulations, the coefficient of restitution for all particle 
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collisions in our inelastic systems, whether with the walls or other hard spheres, 
gives an elasticity of e = 0.99. 

At a lower temperature (T 2 < T c ), particles begin to condense and form a 
solid near the bottom, as shown in Fig. (3). Note that the particles first form 
a loose square-packed lattice, as in Fig. (3), and then progressively evolve into 
a more compact hexagonal lattice structure, as shown in Fig. (4). Once the 
compact hexagonal solid is formed, the density at the bottom remains more or 
less constant at p a — 2^ w .906. This hexagonal structure is then, for the 
most part, permanently retained. For Fig. (3) and (4), the oscillations in the 
solid regime are real, but they are simply the finite size effect, i.e, the hexagonal 
packing in a finite lattice has two more particles in alternative layers. This 
oscillation would disappear in the thermodynamic limit. 

The critical temperature T c is determined by the temperature at which a 
square-packed solid is formed at the bottom layer. Beyond this temperature, 
the solid regime steadily approaches the close-packed hexagonal structure and 
once attained is fixed at p ~ .906. We point out that between T c , where there 
is a square-packed lattice, and lower temperatures, where there is a hexago- 
nally packed lattice, particles squeeze themselves together, expelling holes and 
progressively forming a compact hexagonal solid. Recall that the critical tem- 
perature, Eq.(14), is a function of p and fi . Eq.(14) shows that p itself is 
also a function of po- Therefore, as the system changes from square-packed to 
compact hexagonal packing, so do the values of po, p,Q, and consequently T c . 
This is because as po increases as the solid gets more compact, causing p to 
increase and T c to decrease. 

5 Global Fits to Simulated Data 

Using an Event-Driven Molecular Dynamic simulation written by Stephan 
Luding, we obtained data for height versus density for a two-dimensional vibrat- 
ing bed of hard-discs. Table 1 shows what control parameters we used and how 
we varied them. For each parameter set, we ran several different simulations at 
different temperatures. 



Table 1: Simulation Parameters 



Parameter Set 


Mass (kg) 


Gravity (£) 


Diameter (to) 


M 


1 


8.378 x 10" 9 


9.81 


.0002 


20 


2 


1.676 x 10" 9 


9.81 


.0002 


20 


3 


8.378 x 10" 9 


4.91 


.0002 


20 


4 


8.378 x 10" 9 


9.81 


.0004 


20 


5 


8.378 x 10" 9 


9.81 


.0002 


40 



Using the method of Rcimann Sums in a computer program, we numerically 
calculated the global equation solution from Eq.(5). Then we fit the numerical 
solution for each density profile to obtain the temperature of each system. From 
Fig. (2), (3), and (4), it is clear that the numerical solution to Eq.(5) fits the 



9 



data extremely well. We note that there are slight density oscillations in the 
solid regime for the simulated data due to the presence of the bottom wall. 

Fig. (2), (3) and (4) are all examples of completely elastic systems, with no 
energy dissipation. Fig. (5) and (6) show comparisons between the elastic and 
inelastic trials. For the inelastic case, Fig. (6), it is clear that the Global 
equation still fits the density profile nicely in the liquid-solid regime. However, 
the density profile tends to underestimate the height at lower densities. This 
is because the particles in the fluid regime are active in the sense that they 
are free to move around, undergoing collisions, while those in the solid regime 
are confined in a cage, fluctuating around fixed crystalline positions, and do 
not undergo collisions. This effect is discussed further in [4], where the Global 
Equation was initially developed. 

To calculate the number of frozen layers for each simulation, we use Eq. 
(16), where T is calculated numerically in the simulation and T c is obtained 
using Eq.(14). We find, as expected, that the /i/ decreases with increasing 
temperature. This can be seen in Fig. (7), which is data for the case of initial 
layer number [i = 40, radii r = 0.0001m and mass m — 8.378 x 10~ 9 kg. One 
can see from the linear fit that the data matches the predictions of Eq.(16). The 
theoretical slope calculated from Eq.(16) is —1.631 x 10 11 while the experimental 
slope obtained from the simulated data shown in Fig. (7) is —1.668 x 10 11 . This is 
a percent difference of only 2.45%. The other systems produced similar results, 
all within 5% or less of the predicted values. 

Fig. (8) provides a snap-shot, where T < T c . In this system, we would 
expect to find 30.95 frozen layers. Indeed, we can see from the picture that the 
first 31 layers appear frozen, whereas particles above this layer can move freely. 
These results illustrate validity of the predictions obtained by using the Global 
Equation of State. 

6 Conclusion 

There are points to consider due to the results presented in this paper. First, 
we have demonstrated in this paper, that the Global Equation of State, first 
presented by Luding et al [4] , does indeed account for the liquid-solid transition 
which exists in a system of hard spheres. The profile derived from the equation 
of state is useful for determining the kinetic temperature of the system, from the 
configurational statistics. The breakdown of particle conservation is no longer 
an issue when looking at the density profile, because the Global Equation of 
State, includes those particles frozen at the bottom. However, we can still use 
the concept of a critical temperature to derive the number of frozen layers at the 
bottem, and compare it to what is observed in our simulated data, due to the 
fact that the density becomes constant as the solid regime becomes hexoganally- 
close packed. The region that is close packed is considered the solid regime 
and can be used to determine the number of frozen layers in the system and 
in turn, the thickness of boundary layers. Since only a fraction of grains are 
mobilized under shear [12,13], and avalanches and many interesting dynamics 
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occur in these thin boundary layers [14,15], such a determination should be of 
technological importance. 

Second, the existence of a gravity induced liquid-solid transition of hard 
spheres must have some interesting consequences to higher order kinetic theory, 
in particular with regard to the dynamic behaviors. Unlike particles in the liquid 
regime, those particles in the solid regime are largely confined in cages and fluc- 
tuate around fixed positions. Their motions resemble lattice vibrations rather 
than binary collisions, and it may be a little peculiar, albeit not unphysical, to 
attempt to describe the lattice vibrations using kinetic theory. If so, such a de- 
scription must include much more than binary collisions. However, the Global 
Equation of State gives a complete picture of a system in transition, accounting 
for both the liquid region with binary collisions as well as the solid region with 
the lattice regions. The connection between the two regions, is an empirical one 
derived by Luding et al [4]. This gives us some insight into the kintetic details 
of the transition from liquid to solid regimes, and could be studied further. 

As discussed in the beginning and demonstrated in this paper, this gravity 
induced liquid-solid transition is not a peculiar phenomenon associated with the 
Global Equation of State, but rather an intrinsic transition inherent in a system 
where an excluded volume interaction is dominant. The formation of a solid 
at the bottom is the appearance of a massive occupied low energy state due to 
the Pauli exclusion principle. Therefore, the bridge between the liquid and solid 
regimes, as presented by the Global Equation of State and the upward spread of 
the solid regime should persist because the Pauli exclusion principle is in action 
in real space, even if one may use different approximations [16-25] or may try 
a different form for the pressure, such as the form suggested by Percus-Yevick 
[26], and/or in higher order truncation. It only disappears in the limit when the 
close volume packing density, p becomes one, which is possible only in the case 
of an ideal Appolonian packing [27]. One should also further expand the ideas 
used to derive the Global Equation of State from two to three dimensions, to 
see if a similar equation can be derived and tested. 

Finally, we have shown that the presence of dissipation does not alter the 
condensation picture at all [28,29], as long as the velocity distribution remains 
Gaussian. Previous experiments [30] have demonstrated the non-Gaussian na- 
ture of the velocity distribution, but if the dissipation is small, which is the case 
for the inelastic simulations carried out in this work, the deviation from Gaus- 
sian should be small. We also point out that for hard sphere systems without 
gravity, there exists no typical energy scale, and thus any transition must be en- 
tropy driven, i.e., there exists no critical temperature, and the phase transition 
occurs at a critical volume fraction [31]. However, for the system considered 
in this paper, there exists a typical energy associated with the potential energy 
due to gravity, and thus this transition is not entropy driven, but energy driven, 
and therefore a critical temperature T c must exist. 
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Figure Captions 



Figure 1: This is a system of 1000 elastic two-dimensional spherical particles 
with initial layer number /i = 40, radii r = 0.0001 m and mass m = 8.378 x 10~ 9 
kg. This displays fits to the data using the Maxwell Boltzmann Distribution, 
the Enskog Equation using %4, and Eq.(5), the density profile derived from the 
Global Equation of State. 

Figure 2: This is a system of 1000 elastic two-dimensional spherical particles 
with initial layer number [i — 40, radii r — 0.0001m and mass m = 8.378 x 10 -9 
kg. The temperature of this system is greater than T c , hence all particles are 
fluidized. The red points represents the data points, while the blue is the density 
profile derived from the Global Equation of State. 

Figure 3: This is a system of 1000 elastic two-dimensional spherical particles 
with initial layer number /z = 40, radii r — 0.0001m and mass m = 8.378 x 10~ 9 
kg. The temperature of this system is slightly below T c . The first few layers 
at the bottom are beginning to transition to a solid, while the majority of the 
particles are still fluidized. The red points represents the data points, while the 
blue is the density profile derived from the Global Equation of State. 

Figure 4: This is a system of 1000 elastic two-dimensional spherical particles 
with initial layer number [i — 40, radii r — 0.0001m and mass m = 8.378 x 10~ 9 
kg. The temperature of this system is much less than T c . The large number 
of the bottom layers in the solid-like regime have become hexagonally packed, 
while a smaller number of particles remain fluidized. The red points represents 
the data points, while the blue is the density profile derived from the Global 
Equation of State. 

Figure 5: This is a system of 1000 elastic two-dimensional spherical particles 
with initial layer number /x = 40, radii r — 0.0001m and mass m = 8.378 x 10~ 9 
kg. The global fit for the density profile fits the elastic system of particles 
throughout the entire range of densities. The red points represents the data 
points, while the blue is the density profile derived from the Global Equation of 
State. 

Figure 6: This is a system of 1000 inelastic two-dimensional spherical particles 
with initial layer number \i — 40, radii r — 0.0001m and mass m = 8.378 x 10~ 9 
kg. The global fit for the inelastic density profile works well for high density 
but deviates at lower densities. The red points represents the data points, while 
the blue is the density profile derived from the Global Equation of State. 

Figure 7: This is a graph of layer number versus temperature for a system 
of 1000 clastic two-dimensional spherical particles with initial layer number 
fi = 40, radii r = 0.0001m and mass m = 8.378-E — 9kg. As the temperature 
decreases, more particles become locked in the solid regime. One can see that 
the relationship is linear, as is predicted by Eq.(16). 

Figure 8: This is a snapshot from a system of 1000 elastic two-dimensional 
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spherical particles with initial layer number fi = 40, radii r = 0.0001m and 
mass m = 8.378 x 10~ 9 kg. The temperature of this system is much less than 
T c . This simulation corresponds to the simulated data in Fig. (3). From our 
derivation /if for this system should be 31 layers. Visually, we see that the 
bottom 31 layers, marked in gray, are indeed in a frozen regime, an area where 
particles do not switch position with other particles. 
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